Stationarity

Statistics that allow time series data to be accurately described at all time points require the time series to be stationary. A stationary time series is one whose statistical properties such as mean, variance, autocorrelation, etc. are all constant over time. For a time series to be considered stationary the following conditions must be met:

1. The mean does not depend on time.
2. The variance is finite and does not depend on time.
3. The correlation between data points only depends on how far apart they are in time and not where they are in time.

To estimate means, variance and autocorrelations from a single realization requires us to meet three stationarity conditions. This means the average number of positive flu cases across our realization does not depend on time, or change over time. This would mean that if we repeated the same year with different realizations, those realizations would have the same mean. If we can safely assume a constant mean across all years we can use all the observations to estimate our mean. Similarly if all years have a finite and constant variance across our years we can use all the data to estimate the common variance. The last requirement requires the correlation between data points to be dependent on how far apart they are in time and not where they are in our five year timespan.

From a visual inspection of the realization it is difficult to determine if the mean is constant over time. If the mean were constant over time then every time period would have the same mean, i.e., the mean for December is the same for every realization. Figure 3 plots the mean positive cases by week of the year from 2015 through 2020. The plot shows there are weekly differences in the mean and that the peak for most years from 2015 thorugh 2020 was probably sometime in February. This implies that the expected value of the time series may depend on time and raises an intial concern regarding stationarity, which will be addressed later.

by_week <- most_recent_5_years.flu_positive_cases %>% 
  select(Week, total_flu_cases.positive) %>% 
  group_by(Week)

# get avearge number of cases for the week of each flu season
mean_weekly_cases <- by_week %>% summarise(mean_positive_cases = mean(total_flu_cases.positive))
# plot the mean cases by week for all flu seasons from 2015 through 2020
plot(mean_weekly_cases$Week,mean_weekly_cases$mean_positive_cases, type = "o", 
     main=c(paste("Figure 3"), 
               paste("Mean Flu Cases by Week of the Year from 15 June 2015 through 07 June 2020")),
     cex.main=0.8,
     xlab = "Flu Season Week Number", ylab = "Mean number of Positive Cases")

The correlation between data points (covariance) only depends on how far apart they are in time and not where they are in time (i.e., contstant autocovariance). The ACF plots in Figure 4 and Figure 5 split the data in half to see if the autocorrelations (autocovariance) change over time. Autocorrelations that change over time would imply a non-stationary time series. Comparing the first half of the data to the second half of the data shows the ACFs are not the same. The difference is particularly evident at lag 10. This suggests the autocovariance of the data is not constant over time and a transformation of the data may be required.

# to compare the ACF structure of the first half of the data to the second half.
par(mfrow = c(1,2))
acf(ts.flu[1:130], main="Figure 4 \nAutocorrelations for First Half of Data")
acf(ts.flu[131:260], main="Figure 5 \nAutocorrelations for Second Half of Data")

To deal with the non-constant variance, the data is log-transformed as shown in Fiure 6, which serves to normalize the values.

# add column of log transformed data to datframe
log_data <- most_recent_5_years.flu_positive_cases %>% 
  mutate(log_flu=log(total_flu_cases.positive))
# convert to time series object
ts.log_flu <- ts(log_data[ ,6], start=2015, frequency = 52)

# plot the data
plot(ts.log_flu, main=c(paste("Figure 6"), 
               paste("Log-transformed Weekly Flu Cases in the U.S. from June 15, 2015 through June 07, 2020")), cex.main = 0.9, xlab="Year", ylab="Positive Flu Cases")

A comparison of the ACF plots for the log-transformed data shows the first half (Figure 7) is nearly identical to the second half (Figure 8), which suggests constant autocovariance over time.

# to compare the ACF structure of the first half of the data to the second half.
par(mfrow = c(1,2))
acf(ts.log_flu[1:130], main="Figure 7 \nAutocorrelations for First Half of Log-transformed Data", cex.main=0.8)
acf(ts.log_flu[131:260], main="Figure 8 \nAutocorrelations for Second Half of Log-transformed Data", cex.main=0.8)

ACF and Spectral Density

The ACF plot shown in Figure 9 exhibits converging, sinusoidal behavior, which is characteristic of complex conjugate roots and may suggest an AR(2) process.

# plot the ACF and spectral densities
#invisible allows the plot to print, but supresses the output
acf(ts.log_flu, lag.max = 200, main = "Figure 9")

Additionaly, when looking at the spectral density, which helps identify the frequency content of a time series, there is a significant peak near 0 suggesting complex roots. There appears to be cyclic behavior in the time series data, which could imply seasonality and a non-stationary process. In the first spectral density plot there is a large peak near zero. When the truncation point is changed to 20, this peak is at approimately 0.02, or 1/50, which indicates a period of approximately one year.

# plot the ACF and spectral densities
# invisible allows the plot to print, but supresses the output
par(mfrow = c(1,2))
invisible(parzen.wge(ts.log_flu))
invisible(parzen.wge(ts.log_flu,trunc = 20))

However, it is important to note that time series data with cyclic behavior and no trend or seasonality is considered stationary if the cycles are not of a fixed length.⁴ Cyclic behavior should not be confused with seasonal behavior. If the fluctuations in the data are not of a fixed period, then they are considered cyclic. If the period of the cycle is fixed, then the pattern is seasonal.

From Table 1 below, it is evident that the cycles (measuring from peak to peak) are not of a fixed length and display aperiodic cyclic behavior. The number of weeks that elapse between peaks for the time series varies between 49 and 55 weeks. Intuitively, this makes sense as the peak of the flu “season” doesn’t necessarily fall on the same week or month every year. Therefore, it can be argued that the pattern in the weekly data is not seasonal but cyclic.

Table 1: Flu Season Peak Week

To give some validation to the initial assesments of the data lacking seasonality and trend, a more formal approach to test for stationarity and trend is explored.

Dickey-Fuller Test for Stationarity

Employing an augmented Dickey-Fuller test helps determine if one or more seasonal factors should be included in the model and tests the null hypothesis that the autoregressive model has a root outside of the unit circle. The test depends on failing to reject the null hypothesis to decide whether there is a unit root present. However, it shoul be noted that failing to reject the null hypothesis is not evidence that a unit root (i.e., seasonal factor) exists.

In the case of this data, the augmented Dickey-Fuller test rejects the null hypothesis with a p-value of 0.01, suggesting there are no seasonal factors present and validating the initial inspection of the data.

# Check for stationarity using the Dickey-Fuller test
adf.test(ts.log_flu)  

    Augmented Dickey-Fuller Test

data:  ts.log_flu
Dickey-Fuller = -5.6823, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary

Cochrane-Orcutt Test to Check for Trend

To check if there is a trend, the Cochrane-Orcutt test is employed, which tests the hypothesis that there is no serial correlation in the residuals (i.e., no trend). The Cochrane-Orcutt test yields a p-value of 0.9133 and fails to reject the null hypothesis, which suggests there is no trend.

# Check for trend using Cochrane-Orcutt
df <- log_data[,c(2,6)]
x <- ts.log_flu
t = seq(1,260,1) 
fit = lm(x~t, data = df)

# Cochrane-Orcutt test
#cfit = cochrane.orcutt(fit)
#summary(cfit)

Per the initial visual inspection, the Cochrane-Orcutt test, and the Dickey-Fuller test, the data for the model will be assumed to be stationary with no trend and no seasonal factors.

Model Selection Methodology

The Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are used in this case study to select candidate models. The AIC and the BIC are measures that score a model based on its log-likelihood and complexity. The AIC aims to reduce the white noise variance in the model and penalizes models that add additional terms. The BIC is concurrently employed as the AIC tends to select higher order models, i.e., it may propose selecting an ARMA (2,2) over ARMA (1,1) model. The BIC imposes stronger penalties for increasing the orders of \(p\) and \(q\) and will tend to select models with fewer parameters. Lower values for AIC and BIC are preferred.

The function aic5.wge from the tswge package is used to identify the models with the lowest AIC and BIC. The AIC identifies a potential ARMA(3,1) model.

formattable(invisible(aic5.wge(ts.log_flu, type = 'aic')))

---------WORKING... PLEASE WAIT... 


Five Smallest Values of  aic 

The BIC method also identifies an ARMA(3,1) model as the best structure for the data.

formattable(invisible(aic5.wge(ts.log_flu, type = 'bic')))

---------WORKING... PLEASE WAIT... 


Five Smallest Values of  bic 

Visual Inspection of Residuals

Per the realization plot, the residuals look random, which is characteristic of white noise.

# check residuals for white noise
invisible(plotts.sample.wge(params$res))

The plot of the sample autocorrelations shows lag 25 outside the confidence bands. However, at a 95% confidence level, 1 out of every 20 lags is expected to be outside the bands. Therefore, per a visual inspection, the residuals appear to be white noise.

acf(params$res, lag.max = 60, main="Figure 10")

Ljung-Box Test for residuals

A visual inspection looks at each autocorrelation separately. An alternative approach to checking the residuals for white noise is performing a Ljung-Box test. The Ljung-Box test approaches the autocorrelations as a group to determine if the residuals are white noise. It tests the null hypothesis (\(H_0\)) that all autocorrelations (\(\rho\)) are zero (i.e., the residuals are white noise). \[ H_0: \rho_1 = \rho_2 = ... = \rho_K = 0\]

If at least one autocorrelation is not zero, then white noise is not present. \[H_a: at\;least\;one\; \rho_k \neq 0, \, for\,\, 1 \leq k \leq K\]

In the tswge package, the residuals are found in the output variable $res. These are calculated within the functions est.ar.wge and est.arma.wge.

The Ljung-Box test yeilds \(p > 0.05\), which fails to reject the null hypothesis and suggests the residuals are white noise.

ljung.wge(params$res, p = 3, q = 1)  # pval is > 0.05 and fails to reject the null hypothesis

Obs -0.01099454 0.06186652 -0.1007106 0.03124515 -0.02971618 0.01862058 -0.01032152 0.1088832 -0.04462303 -0.06125367 -0.041286 0.03850463 0.0005598552 0.0218697 -0.002265376 0.03318781 -0.0008419987 -0.0740012 -0.008402974 -0.04907715 0.0265037 0.007055816 0.05382664 0.06023577 

$test
[1] "Ljung-Box test"

$K
[1] 24

$chi.square
[1] 14.7704

$df
[1] 20

$pval
[1] 0.7893909

As a second check, a different K-value is used, which yields \(p > 0.05\) also suggesting the residuals are white noise.

# second check with different K-value
ljung.wge(params$res, p = 3, q = 1, K = 48)  # pval is > 0.05 and fails to reject the null hypothesis

Obs -0.01099454 0.06186652 -0.1007106 0.03124515 -0.02971618 0.01862058 -0.01032152 0.1088832 -0.04462303 -0.06125367 -0.041286 0.03850463 0.0005598552 0.0218697 -0.002265376 0.03318781 -0.0008419987 -0.0740012 -0.008402974 -0.04907715 0.0265037 0.007055816 0.05382664 0.06023577 0.1354727 -0.05804677 0.08891488 0.002882477 0.08403592 0.00791233 -0.06529636 0.02418349 0.007283807 -0.04690817 0.01662693 0.02552653 0.05252992 -0.04462137 -0.03852077 -0.02863463 -0.01039259 -0.01992841 -0.06997014 0.01066953 0.02056869 0.05351268 0.03124399 -0.02713629 

$test
[1] "Ljung-Box test"

$K
[1] 48

$chi.square
[1] 33.38016

$df
[1] 44

$pval
[1] 0.8781805

Per a visual inspection and the Ljung-Box test, the residuals are white noise. This lends confidence in the model selected by the AIC and BIC methods.

Forecasting the Data

A plot of the forecast for the last year of data, 52 weeks, suggests the model does a fairly good job of predicitng the one year’s worth of data. However, this result suggests that the model may be better suited to forecast shorter time horizons.

f = fore.aruma.wge(ts.log_flu,
                   phi=params$phi, 
                   theta=params$theta, 
                   n.ahead = 52, 
                   lastn = TRUE, 
                   plot=TRUE, 
                   limits=TRUE)

Performance Metric

The average squared error (ASE) will be used to measure the goodness of fit of the model (performance). The ASE measure takes the sum of the square of the difference between the predicted value (forecast), \(\hat X_i\), and the actual value, \(X_i\). It then averages the error over \(n\) number of observations. A lower ASE value indicates the model made fewer forecast errors.

\[ASE = \frac{\sum(\hat X_i - X_i)^2}{n}\]

It should be noted that the ASE is a snapshot in time and can vary for the same data set depending on the size of the training data. It uses \(n-k\) values, where \(n\) is the length of the time series, to train the model and then uses the last \(k\) values to validate forecasted values. For the data in this case study, a single ASE for a 6-month (26-week) forecast would fit the model on 235 weeks of data and then test on the last 26 weeks, for example.

A more useful approach is to shorten the training period and fit the model on a smaller training set (a shorter “window” of time) and then validate the data on the subsequent \(k\) values. The training set, or window, then “rolls” or “slides” to the subsequent period (week) and evaluated again and again.

This is called a rolling window ASE. With the rolling window ASE method, we are not simply taking the last ASE observation from \(n\) prior periods. The rolling window characteristic loops through time horizons and averages the ASEs together, which can prove to be a more stable representation of the overall model ASE. For example, if there was some particularly odd behavior in the recent past of a time series, a single ASE could be misleading.

Like cross-validation, we’re taking the model and sliding across the whole dataset and seeing how well it predicts the next \(n\) observations. The training dataset will be comprised of at least 130 weeks, or 2.5 years of data. The forecast horizon will be used as the validation set.

Ideally, for the rolling window ASE charts below, a low and steady ASE value (dotted red line) as compared to the observed values (blue line) is preferred. This would indicate that the model did a good job of predicting most, if not all, the observed values. Spikes in the ASE value represent observed values that were not predicted well, areas of large error.

#Code from Prof. Sadler's Time Series Course Unit 7

#Model 1
phis = params$phi
thetas = params$theta
s  = 0
d  = 0
trainingSize = 130 # this is the window size (we used a window of 2.5 years or 130 weeks)
total_number_of_observations = 260

horizon = 2 # we forecast out 2 weeks
ASEHolder.2_weeks = numeric() # this is an empty varible that will hold all the ASE values

for( i in 1:(260-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(ts.log_flu[i:(i+(trainingSize-1))], 
                             phi = phis, theta = thetas, 
                             s = s, d = d, n.ahead = horizon, plot=FALSE)
  ASE = mean((ts.log_flu[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder.2_weeks[i] = ASE
}

mean_windowed_ASE.2_weeks = mean(ASEHolder.2_weeks)
median_windowed_ASE.2_weeks = median(ASEHolder.2_weeks)

# visualization of windowed ASE over time
par(mar = c(5,5,2,5))
plot(ts.log_flu, type="l", ylab='Flu Positive Rate', xlab='Time', las = 1, col="blue", main = 'Figure 11 \nRolling Window ASE Over Time  (2-week forecast)')
par(new = T)
# plot rolling window ASE 
plot(ASEHolder.2_weeks, type="l", lty=2, axes=F, ylab=NA, xlab=NA, col="red")
# create tick marks and label on right vertical axis
axis(side=4, las=1)
# add ASE line
mtext(side=4, line=3, 'ASE')
# add legend
legend("topleft", legend=c("Obs. Value", "ASE"), lty=c(1, 2), col=c("blue", "red"), cex=.6)

horizon = 4 # we forecast out 1 months, or 4 weeks
ASEHolder.4_weeks = numeric() # this is an empty varible that will hold all the ASE values

for( i in 1:(260-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(ts.log_flu[i:(i+(trainingSize-1))], 
                             phi = phis, theta = thetas, 
                             s = s, d = d, n.ahead = horizon,plot=FALSE)
  ASE = mean((ts.log_flu[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder.4_weeks[i] = ASE
}

mean_windowed_ASE.4_weeks = mean(ASEHolder.4_weeks)
median_windowed_ASE.4_weeks = median(ASEHolder.4_weeks)


# visualization of windowed ASE over time
par(mar = c(5,5,2,5))
plot(ts.log_flu, type="l", ylab='Flu Positive Rate', xlab='Time', las = 1, col="blue", main = 'Figure 12 \nRolling Window ASE Over Time (4-week forecast)')
par(new = T)
# plot rolling window ASE 
plot(ASEHolder.4_weeks, type="l", lty=2, axes=F, ylab=NA, xlab=NA, col="red")
# create tick marks and label on right vertical axis
axis(side=4, las=1)
# add ASE line
mtext(side=4, line=3, 'ASE')
# add legend
legend("topleft", legend=c("Obs. Value", "ASE"), lty=c(1, 2), col=c("blue", "red"), cex=.6)

horizon = 12 # we forecast out 3 months, or 12 weeks
ASEHolder.3_months = numeric() # this is an empty varible that will hold all the ASE values

for( i in 1:(260-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(ts.log_flu[i:(i+(trainingSize-1))], 
                             phi = phis, theta = thetas, 
                             s = s, d = d, n.ahead = horizon,plot=FALSE)
  ASE = mean((ts.log_flu[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder.3_months[i] = ASE
}

mean_windowed_ASE.3_months = mean(ASEHolder.3_months)
median_windowed_ASE.3_months = median(ASEHolder.3_months)

# visualization of windowed ASE over time
par(mar = c(5,5,2,5))
plot(ts.log_flu, type="l", ylab='Flu Positive Rate', xlab='Time', las = 1, col="blue", main = 'Figure 13 \nRolling Window ASE Over Time (3-month forecast)')
par(new = T)
# plot rolling window ASE 
plot(ASEHolder.3_months, type="l", lty=2, axes=F, ylab=NA, xlab=NA, col="red")
# create tick marks and label on right vertical axis
axis(side=4, las=1)
# add ASE line
mtext(side=4, line=3, 'ASE')
# add legend
legend("topleft", legend=c("Obs. Value", "ASE"), lty=c(1, 2), col=c("blue", "red"), cex=.6)

horizon = 26 # we forecast out 6 months, or 26 weeks
ASEHolder.6_months = numeric() # this is an empty varible that will hold all the ASE values

for( i in 1:(260-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(ts.log_flu[i:(i+(trainingSize-1))], 
                             phi = phis, theta = thetas, 
                             s = s, d = d, n.ahead = horizon,plot=FALSE)
  ASE = mean((ts.log_flu[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder.6_months[i] = ASE
}

mean_windowedASE.6_months = mean(ASEHolder.6_months)
median_windowedASE.6_months = median(ASEHolder.6_months)

# visualization of windowed ASE over time
newASE = c(rep(NA, 131), ASEHolder.6_months) # this plots ASE from week 131 onward

par(mar = c(5,5,2,5))
plot(ts.log_flu, type="l", ylab='Flu Positive Rate', xlab='Time', las = 1, col="blue", main = 'Figure 14 \nRolling Window ASE Over Time (6-month forecast)')
par(new = T)
# plot rolling window ASE 
plot(ASEHolder.6_months, type="l", lty=2, axes=F, ylab=NA, xlab=NA, col="red")
# create tick marks and label on right vertical axis
axis(side=4, las=1)
# add ASE line
mtext(side=4, line=3, 'ASE')
# add legend
legend("topleft", legend=c("Obs. Value", "ASE"), lty=c(1, 2), col=c("blue", "red"), cex=.6)

horizon = 52 # we forecast out 1 year, or 52 weeks
ASEHolder.1_year = numeric() # this is an empty varible that will hold all the ASE values

for( i in 1:(260-(trainingSize + horizon) + 1))
{
  forecasts = fore.aruma.wge(ts.log_flu[i:(i+(trainingSize-1))], 
                             phi = phis, theta = thetas, 
                             s = s, d = d, n.ahead = horizon,plot=FALSE)
  ASE = mean((ts.log_flu[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder.1_year[i] = ASE
}

mean_windowedASE.1_year = mean(ASEHolder.1_year)
median_windowedASE.1_year = median(ASEHolder.1_year)

# visualization of windowed ASE over time
par(mar = c(5,5,2,5))
plot(ts.log_flu, type="l", ylab='Flu Positive Rate', xlab='Time', las = 1, col="blue", main = 'Figure 15 \nRolling Window ASE Over Time (1-year forecast)')
par(new = T)
# plot rolling window ASE 
plot(ASEHolder.1_year, type="l", lty=2, axes=F, ylab=NA, xlab=NA, col="red")
# create tick marks and label on right vertical axis
axis(side=4, las=1)
# add ASE line
mtext(side=4, line=3, 'ASE')
# add legend
legend("topleft", legend=c("Obs. Value", "ASE"), lty=c(1, 2), col=c("blue", "red"), cex=.6)

From the rolling window ASE charts above and Table 2 below the performance of the model deteriorates as the horizon increases.

In Table 2, the median and mean rolling window ASE values are displayed

Table 2: Rolling Window ASE

single_ASE.2_weeks = tail(ASEHolder.2_weeks, n=1)
single_ASE.4_weeks = tail(ASEHolder.4_weeks, n=1)
single_ASE.3_months = tail(ASEHolder.3_months, n=1)
single_ASE.6_months = tail(ASEHolder.6_months, n=1)
single_ASE.1_year = tail(ASEHolder.1_year, n=1)

forecast_horizon <- c("2 weeks", "4 weeks", "3 months", "6 months", "1 year")
single_ASE <- c(single_ASE.2_weeks, single_ASE.4_weeks, single_ASE.3_months, single_ASE.6_months, single_ASE.1_year)
mean_rolling_window_ASE<- c(mean_windowed_ASE.2_weeks,mean_windowed_ASE.4_weeks, mean_windowed_ASE.3_months, mean_windowedASE.6_months, mean_windowedASE.1_year)
median_rolling_window_ASE <- c(median_windowed_ASE.2_weeks, median_windowed_ASE.4_weeks, median_windowed_ASE.3_months, median_windowedASE.6_months, median_windowedASE.1_year) 

rolling.window.ASE <- data.frame(forecast_horizon, single_ASE, mean_rolling_window_ASE, median_rolling_window_ASE)

names(rolling.window.ASE)[names(rolling.window.ASE) == "forecast_horizon"] <- "Forecast Horizon"
names(rolling.window.ASE)[names(rolling.window.ASE) == "single_ASE"] <- "Single ASE"
names(rolling.window.ASE)[names(rolling.window.ASE) == "mean_rolling_window_ASE"] <- "Rolling Window ASE (Mean)"
names(rolling.window.ASE)[names(rolling.window.ASE) == "median_rolling_window_ASE"] <- "Rolling Window ASE (Median)"
customGreen0 = "#DeF7E9"
customGreen = "#71CA97"

formattable(rolling.window.ASE, align =c("l","c","c"), list(
  `Forecast Horizon` = formatter("span", style = ~ style(color = "grey",font.weight = "bold")), 
  `Single ASE`= color_tile(customGreen, customGreen0),
  `Rolling Window ASE (Mean)`= color_tile(customGreen, customGreen0),
  `Rolling Window ASE (Median)`= color_tile(customGreen, customGreen0)
))